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ABSTRACT 

We address a primary question regarding the physical mechanism that trig- 
gers the energy release and initiates the onset of eruptions in the magnetar mag- 
netosphere. A self-consistent stationary, axisymmetric model of the magnetar 
magnetosphere is constructed based on a force-free magnetic field configuration 
which contains a helically twisted force-free flux rope. The magnetic field con- 
figurations in the magnetosphere are obtained as solutions of an inhomogeneous 
Grad-Shafranov (GS) equation. Given the complex multipolar magnetic fields 
at the magnetar surface, we also develop a convenient numerical scheme to solve 
the GS equation. Depending on the surface magnetic field polarity, there exist 
two kinds of magnetic field configurations, inverse and normal. For these two 
kinds of configurations, variations of the flux rope equilibrium height in response 
to gradual surface physical processes, such as flux injections and crust motions, 
are carefully examined. We find that equilibrium curves contain two branches, 
one represents a stable equilibrium branch, the other an unstable equilibrium 
branch. As a result, the evolution of the system shows a catastrophic behavior: 
when the magnetar surface magnetic field evolves slowly, the height of flux rope 
would gradually reach a critical value beyond which stable equilibriums can no 
longer be maintained. Subsequently the flux rope would lose equilibrium and the 
gradual quasi-static evolution of the magnetar magnetosphere will be replaced 
by a fast dynamical evolution. In addition to flux injections, the relative mo- 
tion of active regions would give rise to the catastrophic behavior and lead to 
magnetic eruptions as well. We propose that a gradual process could lead to 
a sudden release of magnetosphere energy on a very short dynamical timescale, 
without being initiated by a sudden fracture in the crust of the magnetar. Some 
implications of our model are also discussed. 
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1. INTRODUCTION 

Two intimately connected classes of young neutron stars — Soft Gamma-ray Repeaters 
(SGRs) and Anomalous X-ray Pulsars (AXPs), which are commonly referred to as mag- 
netars, both show high energy emissions (Mazets et al. 1979; Mereghetti & Stella 1995; 
Kouveliotou et al. 1998; Gavriil et al. 2002). It is widely believed that the X-ray lu- 
minosity in these sources is powered by the dissipation of non-potential (current-carrying) 
magnetic fields in the ultra-strongly magnetized magnetosphere with the magnetic field B 
~ 10 14 — 10 15 G (Duncan & Thompson 1992; Thompson & Duncan 1996; Thompson et al. 
2002). Occasionally, a much brighter outburst has been observed, i.e., a giant flare releases 
a total energy of ~ 10 46 ergs and has a peak luminosity of ~ 10 44_46 ergs s" 1 (for recent 
reviews see Woods & Thompson 2006; Mereghetti 2008). Although the energy for magnetar 
outbursts is widely believed to be supplied by the star's magnetic field, the physical process 
by which the energy is stored and released remains one of the great puzzles in high-energy 
astrophysics. Two possibilities exist for the location where the magnetic energy is stored 
prior to an eruption: in the magnetar crust or in the magnetosphere. For the former pos- 
sibility, a giant flare may be caused by a sudden untwisting of the internal (to the neutron 
star) magnetic field (Thompson & Duncan 2001). The subsequent quick and brittle fracture 
of the crust leads to energetic outburst^- During the outbursts, there would be an enhanced 
twist of the magnetospheric magnetic field lines. In this crust scenario, the energy stored 
in the external twist is limited by the tensile strength of the crust. Alternatively, due to 
the difficulties to explain the short timescale of the giant flare rise time, ~ 0.25 ms (Palmer 
et al. 2005), the second possibility — the magnetospheric storage model, was proposed by 
Lyutikov (2006). In this particular scenario, the energy stored in the external twist need not 
be limited by the tensile strength of the crust, but instead by the total external magnetic 
field energy (Yu 2011b). In the magnetospheric storage model, the magnetic energy storage 
processes take place quasi-statically on a longer timescale than the dynamical flare timescale 
prior to the eruption. 

In the magnetospheric model for giant flares, the energy released during an eruption is 
built up gradually in the magnetosphere before the eruption. Some interesting properties 



1 However, recent calculations by Levin & Lyutikov (2012) imply that plastic deformations of the crust are 
more likely to occur and the crust model of giant flares may not explain the fast dynamical energy release. 
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about the storage of magnetic energy of the magnetospheric models have been discussed in 
Yu (2011b). But there still remains a primary question regarding magnetospheric models, 
i.e., what is the mechanism that triggers the energy release and initiates the eruption? More 
specifically, how a very gradual process by the flux injections (Kluzniak & Ruderman 1998; 
Thompson et al. 2002) or crust motions (Ruderman 1991) could lead to a sudden release 
of magnetosphere energy on a very short dynamical timescale, without being initiated by a 
sudden fracture in the rigid component of the neutron star. This catastrophic behavior is 
essentially reminiscent of solar flares and coronal mass ejections (CMEs). It is conceivable 
that the magnetosphere adjusts quasi-statically in response to the slowly- changing boundary 
conditions at the magnetar surface. After reaching a critical point, the magnetosphere could 
no longer maintain a stable equilibrium and a sudden reconfiguration of the magnetic field 
occurs due to loss of equilibrium (Forbes & Isenberg 1991; Isenberg et al. 1993; Forbes & 
Priest 1995). The subsequent physical processes would proceed on a dynamical time scale. 
This catastrophic process naturally explains the puzzle how a very slow process could lead 
to the sudden release of external magnetic energy on a much shorter timescale (Thompson 
et al. 2002). 

The magnetar giant flares may involve a sudden loss of equilibrium in the magneto- 
sphere, in close analogy to solar flares and CMEs (Lyutikov 2006). A number of CMEs show 
structures consistent with the ejection of a magnetic flux ropej, as has been reported by 
Chen et al.(1997) and Dere et al. (1999). Hence magnetic flux ropes have been presumed to 
be typical structures in the solar corona, and their eruptions might be closely related to solar 
flares and CMEs (Forbes & Isenberg 1991; Isenberg et al. 1993). Similarly, in the magnetar 
magnetosphere, magnetic flux ropes could be generated due to the pre-flare activity (Gotz 
et al. 2007; Gill & Heyl 2010). As the magnetic flux injects from deep inside the magnetar, 
the dissipation of the magnetic field may give rise to the precursor activity. The magnetic 
dissipation of the precursor could also lead to topology changes of the magnetic fields and 
the formation of a magnetic flux rop^l Such a flux rope is also an indispensable ingredient 
for the radio afterglow observed in SGR1806 (Gaensler et al. 2005; Lyutikov 2006). It is 
worthwhile to note that the magnetic field interior to the flux rope, which is suspended in 
the magnetosphere, is helically twisted. It corresponds to a locally twisted feature in the 
magnetosphere (Thompson et al 2002; Pavan et al. 2009). Such locally twisted flux ropes 
seem to be more consistent with recent observations, which suggest the presence of localized 



2 The flux rope is a helically twisted magnetic arcade anchored on the solar surface and often used to 
model prominences in the solar corona. 

3 In this work, we do not address the question of how a flux rope might be formed. Possible mechanism 
had been discussed by van Ballegooijen & Martens (1989). 
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twist, rather than global twist (e.g. Woods et al. 2007; Perna & Gotthelf 2008). 

Observations show a striking feature that the emergence of a strong four-peaked pattern 
in the light curve of the 1998 August 27 event from SGR 1900+14, which was shown in data 
from the Ulysses and BeppoSAX gamma-ray detectors (Feroci et al. 2001). These remarkable 
data may imply that the geometry of the magnetic field was quite complicated in regions 
close to the star. It is reasonable to infer that, near the magnetar surface, the magnetic field 
geometry of an SGR/AXP source involves higher multi-poles. The multipolar magnetic field 
configurations could be readily understood within the magnetar model. Physically speaking, 
the electric currents, formed during the birth of magnetar, slowly push out from within 
the magnetar and generate active regions on the magnetar surface. These active regions 
manifest themselves as the multipolar regions on the magnetar surface. Due to the presence 
of the active regions, the magnetic field may deviates from a simple dipole configuration near 
the magnetar surface (Pavan et al. 2009). Our calculations show that multipolar magnetic 
active regions, especially their relative motions, would have important implications for the 
catastrophic eruptions of magnetar giant flares (see Section 4). 

Motivated by the similarity between giant flares and solar CMEs, Lyutikov (2006) spec- 
ulated that magnetar giant flares may also be trigged by the loss of equilibrium of a magnetic 
field containing a twisted flux rope. But no solid calculations about the equilibrium loss of a 
flux rope in magnetar magnetosphere have been performed yet. In this paper we focus on the 
possibility of magnetospheric origin for giant flares and propose that the gradual variations at 
the magnetar surface could lead to fast dynamical processes in the magnetosphere. We will 
construct a force-free magnetosphere model with a flux rope suspended in the magnetosphere 
and study the catastrophic behavior of the flux rope in a background multi-polar magnetic 
field configuration, taking into account the possible effects of flux injections (Kluzniak & 
Ruderman 1998; Thompson et al. 2002) and crust horizontal motions (Ruderman 1991). 
We are especially interested in the critical height of the flux rope that can be achieved in 
our model. In the mean time, we also develop a convenient numerical scheme to solve the 
inhomogeneous Grad-Shafranov (GS) equation. Since observed magnetars have a very slow 
rotation rate, we ignore rotation effects throughout this work. 

This paper is structured as follows: in Section 2 we describe the basic equations for 
the force-free magnetosphere model as well as the multipolar boundary conditions. Two 
possible magnetic configurations are also discussed in this section. In Section 3 we will 
discuss the internal and external equilibrium constraints in our model. Numerical results 
about catastrophic behaviors of the magnetosphere in response to flux injections and crust 
motions are discussed in Section 4. Conclusions and discussions are given in Section 5. 
Technical details about the force-free magnetosphere magnetic field are given in Appendix 
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A and B. 



2. Axisymmetric Force-Free Magnetosphere with Multipolar Boundary 

Conditions 

The magnetic fields in the magnetar magnetosphere are so strong that the inertia and 
pressure of the plasma could be ignored (Thompson et al. 2002; Yu 2011a). As a result, the 
magnetosphere is assumed to be in a force-free equilibrium state, in which J x B = 0. The 
axisymmetric force-free magnetic field configurations is determined by an inhomogeneous 
Grad-Shafranov (GS) equation. Throughout this paper, we work in the spherical polar 
coordinates (r,0,(j>). 

2.1. Force-Free Magnetic Field Containing a Flux Rope 

In our magnetosphere model, one of the distinguishing features is that there exists a 
helically twisted flux rope in the magnetosphere. The precursor of a giant flare could be 
relevant to the formation of such helically twisted flux ropes (Gotz et al. 2007; Gill & Heyl 
2010). Due to the presence of the flux rope, the magnetic fields consist of two parts, one is 
the fields that are inside the flux rope and the other is the fields outside the flux rope. 

In our model the magnetic twist of the flux rope is locally confined to the flux rope 
interior. This is quite different from the globally twisted magnetic field configurations in 
Thompson et al. (2002) and Beloborodov (2009). They considered a non-potential force- 
free field where the electric currents permeate through the entire magnetosphere, while our 
model only contains a electric current in a spatially confined region, i.e., interior to the flux 
rope. Note that the toroidal flux rope has a minor radius, r , which is small compared 
to the height of flux rope, h, which is actually the major radius of the flux rope. Under 
such circumstances, the magnetic field produced by the current inside the flux rope can be 
viewed as that produced by a wire carrying the net current / at the center of the flux rope 
(Forbes & Priest 1991) and a simple Lundquist (1950) force-free solution could be applied 
to represent the distribution of current density and magnetic field inside the flux rope. A 
brief yet self-contained description of the Lundquist solution is given in Appendix A. 

Outside the flux rope, the magnetic field is essentially potential, i.e., the field outside 
the flux rope is non-twisting. In the regions exterior to the flux rope, the steady state 
axisymmetric magnetic field in the magnetosphere has only poloidal components and can be 
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written as 

B = WxV0, (1) 

where \l/(r, 6) is the magnetic stream function and (f> is the third component of the spherical 
polar coordinates. Written explicitly, the magnetic field is 

B= -1_ (l d 4, ™). (2) 



rsin# \r 86 dr J 

The force- free condition can be cast into the standard Grad-Shafranov (GS) equation (Thomp- 
son et al. 2002) 

d 2 ^ sin6 d ( 1 , n Aix 

+ -^~EE —n^ = -(rsme)-J^ , (3) 



where c is the speed of light. The current density on the right hand side of the above 
equation is caused by a toroidal force-free magnetic flux rope mentioned above, which is 
suspended in the magnetosphere at the height, h, by force balances. Note that the stream 
function ^ is determined simultaneously by the electric current inside the flux rope and 
boundary conditions at the magnetar surface. The boundary conditions will be discussed 
separately in next section. The electric current inside the flux rope could be treated as a 
source term on the right hand side of the above inhomogeneous GS equation. We treat the 
current density of the flux rope as a circular ring current of the form (Priest & Forbes 2000) 

J^= I -5( C os6)5(r-h) , (4) 
h 

where / is the electric current carried by the flux rope. Similar treatments have been adopted 
in the coronal mass ejection (CME) studies (Forbes et al. 1991; Lin et al. 1998). It is clear 
from this equation that the flux rope is located at the equatorial plane {6 = tt/2) and the 
flux rope is the only current source in the region r > r s , where r s is the magnetar radius 
(also see Fig. 2). 



2.2. Multipolar Boundary Conditions at Magnetar Surface 

In order to solve the boundary-value problem associated with the inhomogeneous GS 
equation (J3]), we still need to know the boundary condition at the magnetar surface r = r s 
(where r s is the magnetar radiusjE We choose \l/ at the magnetar surface r = r s to be 

* s (r s ,/i) = ^ aQ(fi) , (5) 



The boundary condition at r — > oo is simply |V$| — > 0, which is satisfied trivially in this work, see 
Appendix B. 
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where the subscript s denotes the value on the neutron star surface, \l/o is a constant with 
magnetic flux dimension, a is a dimensionless quantity that determines the magnitude of 
the flux at the surface, and fi = cos 6. The magnetic field configuration of neutron stars 
is basically a dipole field. But near the neutron star surface, where the loss of equilibrium 
occurs, the magnetic field presents much more complex behaviors (Feroci et al. 2001). To 
simulate multipolar regions on the neutron star surface, we add two Gaussian functions to 
the usual dipole field and consequently the function 0(/i) in the above equation takes the 
following form, 



90) = (1- fi 2 )+exp 



2w 2 



+ exp 



(fjL + IM))' 

2w 2 



(6) 



where /io and w are parameters that determine the magnetic flux distributions at the neutron 
star surface. We take w = 0.001 throughout this paper. It is worthwhile to note that, 
according to the parameter \l/o introduced in Equation (jSJ), a dimensional current can be 
defined as 

4 = ^, (?) 

r s 

where r s is the radius of the magnetar. Throughout this paper, we scale all lengths by the 
neutron star radius r s , magnetic flux by \l/o and current by Io = ^oc/r s . We also define a 
dimensionless current J = I/I for later use, where I is the electric current carried by the 
flux rope (see Section 2.1). 

Fig. [TJ shows the profile of the flux function 0(/i) in Equation ([6]). Note that the 
derivative with respect to /i gives the radial component of magnetic field at the magnetar 
surface. This boundary flux distribution is symmetric with respect to the equator 9 = it/ 2 
or fi — 0. In real circumstances, the active regions on magnet ars may form much more 
complicated patterns without any symmetry. For simplicity, we focus in this paper only on 
systems with such symmetry. The distance between the two active regions is specified by 
the parameter /io- The distribution of 0(/i) for three different values of /io = 0.1,0.2, and 
0.3 are shown in three panels of Fig. 1. It is clear from this figure that, with the increase 
of /io, the active regions move away from each other. Another important physical process, 
flux injections, can be interpreted in terms of the variation of the parameter, a, in Equation 
(jSJ). This kind of variations do not change the shape of 0(/i), but change the magnitude 
of flux. For instance, if an opposite polarity magnetic flux is injected from below, due to 
the magnetic cancellation with the pre-existed magnetic flux at the magnetar surface, the 
absolute value of the parameter |a| might decrease. Some interesting consequences from 
both kinds of alterations will be explored in this paper (see Section 4). 
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Fig. 1. — Distribution of magnetic flux at the magnetar surface. Three panels correspond 
to Hq = 0.1, fiQ = 0.2, and /i = 0.3, respectively. The increase of /j,q clearly shows that the 
magnetic active regions are moving apart, which will have significant implications for the 
catastrophic behavior of flux rope. 
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2.3. Inverse and Normal Configurations 

The solution to Equation (J3J) associated with the boundary condition (jSJ) is of vital 
importance for our further discussion. With some trivial boundary conditions, this equation 
can be solved analytically by the Green-function method (Lin et al. 1998). When more 
complex multipolar boundary conditions, such as Equations ([S]) and (ED, are introduced, this 
GS equation can generally be solved by the variable separation method. In this paper, we 
develop a numerical method to solve this equation and full solutions to this inhomogeneous 
GS equation are given in details in Appendix B. 

When solving the GS Equation (j3J) together with the boundary condition Equations (J5J) 
and OH]), we find that there exist two kinds of magnetic field configurations in the magne- 
tosphere. One is the the normal configuration, which means that the magnetic field at the 
position (r, 6) = (h — ro, 7r/2) threads across the flux rope in the same direction as the 
magnetar surface magnetic field underneath at the equator (r, 9) = (r s , tt/2). The other is 
the inverse configuration, in the sense that the magnetic field at the position (r, 9) = (h — r , 
7r/2) threads across the flux rope in the opposite direction to the magnetar surface magnetic 
field underneath at the equator (r, 9) = (r s , tt/2). In our calculations, the current J is al- 
ways kept positivd^l So for our particular choice of boundary conditions, if a is negative, we 
will get an normal configuration, if a is positive, an inverse one is obtained. Two schematic 
figures, both inverse and normal, are shown in the left and right panel of Fig. 2, respectively. 

3. Local-Internal and Global- External Equilibrium Constraints 

In what follows, we consider that the magnetar magnetosphere evolves on a sufficiently 
long timescale so that we can treat the magnetosphere as being essentially in a quasi-static 
equilibrium. The condition for the flux rope equilibrium includes two parts: the local-internal 
equilibrium and global-external equilibrium (Forbes & Isenberg 1991). 

3.1. Local-Internal Equilibrium Constraints 

For the local internal equilibrium, we assume that the force- free condition, J x B = 0, 
also applies within the flux rope. We adopt the Lundquist (1950) force- free solution to 



5 This is to avoid the negative values of flux rope minor radius ro, according to Equation ((5J. 



-2.5 



-2.5 







0.5 



1.5 



2 



2.5 







0.5 



1.5 



2 



2.5 



x/r 



x/r 



s 



s 



Fig. 2. — Magnetic field lines for both the inverse (left panel) and normal (right panel) 
magnetic field configurations. Field lines are obtained with solutions of the inhomogeneous 
GS equation (see Appendix B). The thick black line denotes the magnetar surface. The 
dashed line is a circle with a radius r = h, where h is the height of the flux rope. The flux 
rope lies at the position (r, 6) = (h, n/2). At the magnetar surface r = r s , two additional 
active regions appear due to our choice of boundary conditions. 
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represent the distribution of current density and field inside the flux ropqj. For a circular 
toroidal flux rope in our case, the Lundquist solution inside the flux rope is still valid as long 
as the minor radius r is much smaller than the major radius, h (also know as the flux rope 
height). In this case, a simple relation between the minor radius of the flux rope r and the 
current flowing in the flux rope, J, can be written as 

r ooA) r 00 /q n 

r ° = — = T' (8) 

where J is the dimensionless current scaled by Jo and r 0Q is the value for r as J = 1. We 
take roo = 0.01 throughout this paper. It should be noted that the internal equilibrium 
constraint, Equation (|E|), is actually an alternative way indicating the conservation of axial 
magnetic flux inside the flux rope (see Appendix A). 



3.2. Global-external Equilibrium Constraints 

The global equilibrium is satisfied when the total force exerted on the flux rope vanishes. 
The ring current inside the flux rope provides an outward force. Intuitively, the anti-parallel 
orientation of the current flowing on the opposite sides of the ring produces a repulsive 
force similar to the force between two parallel wires with anti-parallel currents does. The 
magnitude of this force is equal to the current, J, times the magnetic field B s (Shafranov 
1966): 

I (. 8h 



where r is the minor radius of the toroidal flux rope. The additional terms in the bracket 
of the above equation appear due to the curvature effects of the circular ring current 0. This 
ring current induced force must be balanced by the external field B e . The external magnetic 
field B c at r = h and 6 = rr/2 can be written explicitly as (The contribution from the current 
inside the flux loop must be excluded for the external magnetic field B e . For details, see 
Appendix B) 

B C =J2 nY n d n h~ n - 2 , (10) 

n odd 



6 Thc Lundquist solution is obtained in cylindrical coordinates (see Appendix A). Strictly speaking, the 
Lundquist solution is only applicable for a straight cylindrical twisted flux rope. 

7 For two straight wires, terms in the bracket disappear and the induced magnetic field is strictly propor- 
tional to the electric current and inversely proportional to the distance between the wires. 
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where the coefficients r„'s and d n 's are all given explicitly in Appendix B. The mechanical 
equilibrium condition, which requires matching the external field B e with B s , reads 



where 



f(a,J,h) 



f(a,J,h)=0 



nY n d n h- n - 

,n odd 



J 
h 



TOO 



(11) 



(12) 



Note that the local-internal equilibrium constraint has also been exploited in this equation. 
Note that the scaled current J = I/I Q (not the current J) appears in this equation. This is 
because all quantities, such as h, d n , and J, in this equation are measured in the dimensional 
units mentioned above (see Section 2.2 and also Appendix B). 

The stream function \1/ satisfies the ideal frozen-flux condition, which provides a link 
between the electric current flowing inside the flux rope and the boundary conditions at the 
neutron star surface. Specifically, it requires that the stream function on the edge of the flux 
rope remain constant as the system evolves. At the equator 9 = ir/2, the edge of the flux 
rope is located at r = h — tq and the frozen-flux condition can be written explicitly as 



7T 



const 



(13) 



where h and r are the major and minor radius of the flux rope, respectively. Substituting 
r = h — ro and = n/2 into the stream function, i.e., Equation (Bl) in Appendix B, we 
arrive at another constraint, 

g(a, J, h) = const , (14) 



where 



g(a,J,h) = 



n odd 



r„c n ( i 



Jh) 



n+l 



+ T n d n { h 



rboY 
J J 



(15) 



where d n 's and T n 's have the same meaning as those in Equation (fT2l) . c n 's are explicitly 
given in Appendix A. Note again that all quantities in this equation are measured in the 
dimensional units defined above. The quantity r 00 appears because of the internal equilibrium 
constraint. 

In summary the equilibrium constraints, including the force balance and the frozen-flux 
condition, can be written as the following form 



f{a,J,h) = 
g(a, J, h) = const 



(16) 



where functions / and g are defined in Equations (j!2p and ffT5l) . Numerical values of / and 
g are calculated following the procedures presented in Appendix B. For a given value of a, 
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the above two equations could be treated as a nonlinear set of equations for J and h, which 
can be solved numerically by the Newton- Raphson method (Press et al. 1992). Note that 
the Jacobi matrix necessary for the Newton-Raphson method is hard to obtain analytically 
and we calculate the Jacobi matrix numerically instead. 

4. Loss of Equilibrium in Response to Variations at Magnetar Surface 

We consider the possibility that the primary mechanism for driving a magnetar giant 
flare is a catastrophic loss of equilibrium. The loss of equilibrium is initiated by slow changes 
at the magnetar surface. Physically, there are generally two possible processes that could 
occur at the magnetar surface. One is that new magnetic fluxes, driven by the plastic defor- 
mation of neutron star crust, may be injected continuously into the magnetosphere (Kluzniak 
& Ruderman 1998; Thompson et al. 2002; Lyutikov 2006; Gotz et al. 2007). Another in- 
teresting possibility is brought about by the crust horizontal movement (Ruderman 1991; 
Jones 2003). It is very difficult to compress magnetar crust material very much, or to move 
elements of crust up or down. It is, however, much easier to move parts of the crust hor- 
izontally, in ways which apply only shear strains to it (Thompson & Duncan 2001; Jones 
2003). It is possible that, when the magnetic field is strong enough, the interior magnetic 
stress may cause the active regions of the crust to move horizontally (Ruderman 1991). 

For the first possibility, as the new current- carrying magnetic fluxes are injected, a 
direct consequence is that the background magnetic field would vary gradually because of the 
active flux injections prior to large outbursts. The background magnetic field would increase 
(decrease) if the same (opposite) polarity flux is injected. Variations of the equilibrium 
height of flux rope with alterations in the background magnetic field are carefully examined 
for both the inverse and normal magnetic configurations. In this case, we fix the value of 
fio — 0.1 and investigate the effects of variations of a on the flux rope equilibrium height. 
Numerical results of Equation f fl6|) are shown in Fig 3. and Fig 4. These two figures show 
the results for the normal and inverse magnetic configurations, respectively. The curves 
in Fig. 3 and Fig. 4 consist of two branches, which comes from the fact that there exist 
two roots of h for each particular value of a and the two roots lie on separate branches. 
The upper branch denotes an unstable equilibrium state because when the equilibrium is on 
the upper branch, an slight upward vertical displacement will generate an outward driving 
force. The lower branch, however, stands for a stable branch, in the sense that an slight 
upward displacement would create an inward restoring force, just like an harmonic oscillator. 
The stability of the flux rope could be understood in terms of a spring model. The spring 
coefficient of Hooke's law determines the stability of the spring. It would be instructive to 
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treat the total force, T, as a function of flux rope height, h, while keeping the flux frozen 
condition satisfied. Detailed analysis shows that derivative dT/dh (which is equivalent to 
the spring coefficient of Hooke's law) is negative if the flux rope lies on the lower branch and 
positive on the upper branch (see Fig. 6.18 in Forbes 2010). Negative dT/dh corresponds to 
a normal Hooke spring coefficient and a stable equilibrium, while positive dT/dh corresponds 
to an anomalous Hooke spring coefficient and an unstable equilibrium. The two branches 
are connected by a critical point (nose point). The instability threshold lies at the nose 
point. The nose point can also be understood as the critical loss-of-equilibrium point (red 
point in Fig. 3 and Fig. 4). Once the equilibrium reaches the loss-of-equilibrium critical 
height, the system would no longer stay in a stable equilibrium state. The flux rope will 
lose equilibrium and lead to an eruption. In the case of normal configuration, Fig. 3 shows 
that the increase of the parameter \a\ (flux injection of the same polarity) would bring the 
system to the critical point, which means that the enhancement of the background flux would 
trigger the catastrophic behavior for a normal configuration. While for the case of inverse 
configuration, Fig. 4 shows the decrease of the parameter |er| (flux injection of the opposite 
polarity) would lead to loss of equilibrium, indicating the decay of the background flux works 
for an inverse configuration. Our calculations show that the critical height for the two kind 
of configurations differs much. The normal configuration shows a rather low critical height, 
roughly 3% percent above the magnetar surface, which is, for a typical neutron with radius 
10 6 cm, roughly 3 x 10 4 cm. While for the inverse configuration, the critical height is about 
20% above the magnetar surface, which is about 2 x 10 5 cm. Given the regular arrangements 
that occur at the magnetar surface, the small critical height of the normal configuration 
would indicate that it may not survive those arrangement at the magnetar surface and the 
inverse configuration, whose critical height is larger, is preferred in real circumstances. 

As the interior magnetic stress by the ultra-strong magnetic field in magnetar may 
cause active regions of the crust to move horizontally (Ruderman 1991; Jones 2003; Lyutikov 
2006), the relative positions between multipolar active regions may vary, moving apart or 
approaching each other, which may also have significant implications for the catastrophic 
behavior of magnetospheres. Consequently, we further investigate the response of the flux 
rope to horizontal motions of active regions at the magnetar surface. In our simplified 
model, the distance between two active regions is determined by a single parameter fj, Q in 
Equation fl6]). The increase of /io may indicate that active regions move apart (see Fig.l). 
To investigate the effects of horizontal motions, we fix the value of a and vary the parameter 
/io. Again, we numerically solve Equation (TIB"]) and get two roots of h for each particular 
value of /io- We show in Fig. 5, taking the inverse configuration as an example, the variation 
of the equilibrium height with the distance between two active regions. Similar to Fig. 3 
and Fig. 4, the upper branch and the lower branch denotes unstable and stable equilibrium, 
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Fig. 3. — Equilibrium height h as a function of a . This curve is numerically obtained 
as solutions to Equation ( ITSj) . We always choose the current to be positive. The negative 
value of o denotes a normal magnetic configuration. In this case, with the increase of of |cr|, 
the equilibrium height gradually increases and reaches the critical loss-of-equilibrium point, 
beyond which the flux rope could not maintain the stable equilibrium. The critical height 
for the normal magnetic configuration is approximately h c = 1.028. All lengths are scaled 
by ?v 
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Fig. 4.— The same as Fig. 3 but for an inverse magnetic field configuration. In this 
case, with the decrease of of \a\, the equilibrium height gradually increases and reaches 
the critical loss-of-equilibrium point, beyond which the flux rope could not maintain the 
stable equilibrium. The critical height for normal magnetic configuration is approximately 
h c = 1.24. All lengths are scaled by r s . 
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respectively. As the two active regions move apart, the equilibrium height gradually increases 
and reaches the critical height h c = 1.32 when /x = 0.14. After reaching the critical point, 
the system would no longer maintain a stable equilibrium state. This means that, in addition 
to flux injections, the horizontal motions of the active regions could give rise to the loss of 
equilibrium and dynamical eruptions as well. 

5. Conclusions and Discussions 

In this work, we consider the possibility that the primary mechanism for driving an 
eruption in magnetar giant flare is a catastrophic loss of equilibrium of a helically twisted 
flux rope in the magnetar magnetosphere. The loss of equilibrium behavior of a flux rope is 
investigated in a multi-polar magnetic field configuration, taking into account possible effects 
of flux injections and crust horizontal motions. The loss of equilibrium model describes a 
quasi-static equilibrium that varies in response to slow changes at the magnetar surface. 
Beyond a critical point, the stable equilibrium can not be maintained and the transition to 
a dynamical evolution naturally occurs. 

Equilibrium states of a stationary, axisymmetric magnetic field in the non-rotating mag- 
netosphere containing a flux rope are obtained as solutions of the inhomogeneous GS equation 
in a spherical polar coordinate. In view of the complex multipolar boundary conditions at 
the magnetar surface, we develop a numerical method to solve the GS equation. Two kinds 
of magnetic field configurations, inverse and normal, are carefully examined in this work. 
Both of them present the loss of equilibrium behavior. We carefully examined the critical 
height of the flux rope beyond which a stable equilibrium could not be maintained and a 
sudden release of magnetosphere will be triggered. We find that the critical flux rope height 
is different for the two types of configurations. We also investigate effects of another form of 
boundary changes, crust horizontal motions, on the loss of equilibrium behavior. Our results 
show that both the flux injection and crust motions could trigger the catastrophic behavior 
in the magnetosphere. 

In our simplified model, the flux rope is assumed to be a closed current ring encircling 
the magnetar. It is suspended in the magnetar magnetosphere and two ends of the flux rope 
are not anchored to the magnetar surfaces. We expect that the overall catastrophic behavior 
of our model should remain the same even the anchoring effects of flux rope is taken into 
account. In order to further understand the anchoring effects on the catastrophic behavior, 
more realistic three dimensional model that includes a flux rope with two ends anchored to 
the magnetar surface is worth further investigation. 
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Fig. 5. — Equilibrium height h as a function of fiQ. We fix a = 8.5 in this figure. With 
the increase of of /xq, the two active regions move apart. The equilibrium height gradually 
increases and reaches the critical point, beyond which no stable equilibrium state exists. The 
critical height for normal magnetic configuration is approximately h c = 1.32. All lengths are 
scaled by r s . 
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The magnetic energy that could be released in our model is an interesting issue that 
is worthwhile to be explored further. For a spherical coordinate in our paper, the ground 
energy reference state, based on which the fraction of released magnetic energy can be 
calculated, is the fully open Aly-Sturrock field (Aly 1984, 1991; Sturrock 1991; Yu 2011b). 
Since the boundary conditions in our paper is more complex (not dipolar or quadrupolar 
boundary conditions), the construction of the Aly-Sturrock field is technically nontrivial. 
As a result, the fraction of the energy that could be released in our model involves a very 
careful calculation of the Aly-Sturrock field. Intuitively, according to prior studies, which 
have shown that magnetic configurations with more complex boundary conditions would be 
able to release more energies (Forbes & Isenberg 1991; Isenberg et al. 1993), as well as the 
complex boundary conditions adopted in our model, it is expected that our model would 
release enough energy for a magnet giaut darfl Full details of the energetics of our model 
would be discussed elsewhere. 

It is possible that the current sheet forms after the system loses equilibrium (Forbes 
& Isenberg 1991; Forbes & Priest 1995). With the formation of current sheet, the tearing 
instability would develop inside the current sheet (Komissarov et al. 2007) and the subse- 
quent magnetic reconnection would further accelerate the flux rope (Priest & Forbes 2000). 
Magnetic field configurations with the current sheet in a spherical polar coordinate is a long- 
standing unresolved problem. We note that numerical method developed in this work can 
be further extended to allow the presence of current sheets (Yu 2011b). Further discussions 
about the current sheets formation and their effects on the catastrophic behavior would be 
reported in a separate paper. 

It would be interesting to to examine the spectral properties of the model presented in 
this paper (Thompson et al 2002), in which a locally twisted flux rope is self-consistently 
incorporated into the magnetar magnetosphere. By fitting the spectral features with obser- 
vations (Pavan et al. 2009), certain parameters in this model, e.g., flux rope height, electric 
current and magnetic field, may be better constrained. In parallel, recent Fermi observation 
of Crab nebulae gamma-ray flare could possibly be explained by the magnetic reconnection 
models (Abdo et al. 2011), in which the loss of equilibrium may be the trigger for the 
formation of current sheet and subsequent magnetic reconnection processes. The high en- 
ergy flare emission from Crab Nebulae is thought to be synchrotron radiation by relativistic 
electron-positron pairs accelerated in this current sheet (Uzdensky et al. 2011). It would be 



8 Typically 1% of the magnetic energy release could already account for a giant flare. For a simple dipolar 
boundary condition, about 1% of the magnetic energy can be released (Forbes & Isenberg 1991). However, 
about 5% of the magnetic energy can be released for a quadrupolar boundary condition (Isenberg et al. 
1993). 
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instructive to calculate, based on the model shown in this paper, the emission spectra and 
compare them with recent Crab Nebula observations. 

The models constructed in this work are likely to be useful as initial states in high reso- 
lution force-free electrodynamic numerical simulations to explore the dynamics of magnetic 
eruptions (Yu 2011a). Our current model can not address the dissipation processes that 
occur during giant flares. This is left for a future work to directly simulate the behavior 
of loss of equilibrium and relevant dissipation processes using a newly developed resistive 
force-free electrodynamic code (Yu 2011a). 

Discussions with T. Forbes & J. Lin are highly appreciated. We are grateful to the 
anonymous referee's insightful comments that improve this paper. The research is supported 
by the Natural Science Foundation of China (grants 10703012, 11173057), the Western Light 
Young Scholar Program. The computation in this work is performed at HPC Center, Kun- 
ming Institute of Botany, CAS, China. 
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A. Lunquist's solution Interior to Flux Ropes 

When the flux rope minor radius is small compared to the major radius of the flux rope, 
the flux rope interior magnetic field could be approximated as a straight force-free magnetic 
cylinder. We use Lundquist's solution to represent the flux rope interior magnetic field. The 
solution satisfies the force-free condition V x B = A B, where A is a constant. The solution 
can be written in a cylindrical coordinate (r, 9, z). Note that the z direction in this Appendix 
is actually the <fi direction in the main text. Explicitly, Lundquist's solution is 

B r = , (Al) 

Be^BoJ^Xr), (A2) 

B z = B J {Xr) , (A3) 

where B z is along the central axis, Bq is the azimuthal component, and B r is the radial 
component, Jo and J\ are the zeroth- and first-order Bessel functions. The conservation of 
the toroidal flux simply gives 

"2tt pro 

B z rdrd9 = const . (A4) 



'o jo 

Substituting the force-free condition \B Z = 4iij z /c, we have that 

4tt f 27T r ro , ,„ 47r/ 



, . j z rdrd9 = — — = const , (A5) 
Ac J J Ac 

which can be written equivalently as 

— = const . (A6) 
A 

At the surface of the flux rope, B z is zero so Jo(Ar ) = 0. Therefore Ar is the first zero of 
J and Ar = 2.405. We can finally arrive at 

r I = const , (A7) 

which is Equation (jSJ) in the main text. 



B. General Solution of the inhomogeneous Grad-Shafranov Equation with 

Multipolar Boundary Conditions 

According to the variable separation method, the general solution to the Grad-Shafranov 
equation can be written as (see Yu 2011b) 

V=Y1 (cnRn(r) + d n r~ n ) 

n odd 



n-1 



n+1 



2n + 1 



odd 



(Bl) 
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where c n and c? n are constant coefficients to be specified, P n _i(/i) and P n+ i(/j,) are Legendre 
polynomials, = COS6 1 . Here the piece-wise continuous function R n (r) in the above equation 
is defined as 

{r/h) n+1 r<h 
(h/r) n r > h 

Note that the derivative of the function R n (r) is discontinuous. This feature is exploited in 
the following to handle the inhomogeneous source terms associated with the Dirac-5-type 
current density. Also note an identity for the Legendre polynomials 



Rn(r) 



(B2) 



P n ^(l2) - P n+1 (fl) = (l-/i 2 ) dP n 

2n + 1 n(n + 1) d\x 

The inhomogeneous Grad-Shafranov (GS) equation reads 

1 0^ 



2 ^ sin9 

Or 2 



r 2 06 \ sin 6 06 



r sin 6 — Js . 

c 



Substituting Equation (Bl) and the explicit expression of (Equation 
mogeneous GS equation, we arrive at 



(B3) 
(B4) 

into the inho- 



(1 - /i 2 ) dp n 

n(n + 1) d\i 



d 2 R n 
dr 2 



n(n + 1 — 



-r sin 6 



AtcI 
he 



5{cos6)5{r-h) . (B5) 



Integrating r over an infinitesimally thin shell around r = h, we can rewrite the above 
equation as § 



{i-^) dp n 

n n{n + 1) d\i 



■(2n + l)i 



AttI 



sin 68 (cos 1 



(B6) 



Multiplying sin 6^- on both sides of the above equation and integrating 6 over [0,7r] , we 
have that 



;-l)Vn! 



4nlh 



_2»(==i)!(2=i)!_ 

It is more convenient to calculate c n numerically by the following recursive relation, 

'n + 2 s 



fAnIh\ 
ci = 0.5 I — - — 1 , c n+2 



n + 1 



(B7) 



(B8) 



Note that terms in the stream function involving c„'s are induced by the current inside the 
flux rope and the external magnetic field in Equation (I1UI) only involves terms related to d n 's 



9 Note the fact that the first order derivative of R n (r) is discontinuous, and ^r\ r _ h+ — ^r\ r -h- = 
-(2n + l)i. " 
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in Equation (Bl). According to Equation (j2J), the external magnetic field B c at r = h and 
= 7r/2, can be written explicitly as 

5e= E nr n rf n /l- n - 2 , (B9) 



n odd 



where 



P n -i(M = 0)-P„+i(M = 0) 



2n + 1 

The coefficient T n can be readily calculated as 

Ti = 0.5 , T n+2 = — 



n 



n 



(BIO) 



(Bll) 



To specify the coefficients of d n 's, we have to take into account the boundary conditions 
at the magnetar surface. We require that the stream function \I/ be equal to i.e., the 
boundary conditions (jHJ) at the neutron star surface r = r s . To achieve this, we expand 
\l/ s (/i) as the following 

"P n _i(//) - P n +iO)' 



The coefficients a n 's can be expressed as 

2n + 1 



2n + 1 



1 0(0) . 
i an 



Finally the coefficients d n 's can be obtained as follows, 

oL = r! 



c„. — 



T \ n+1 
1 s i 



(B12) 



For numerical conveniences, we scale all lengths by r s , the magnetic flux by \l/o and the 
current by Jo = ^oc/r s throughout this paper. Then the equation B e = B s becomes 

f(a,J,h)=0, (B13) 

where 



-n-2 



J 



h 



^00 



/(a, J,h)= £ nT n d n h 

n odd 

and J = I/Iq . Similarly the frozen-flux condition can be written as 
g(<r, J,h)= £ T n c n (l - j£j 

n odd 



r„(/„ ( h - 1 -y 



const . 



(B14) 



(B15) 



Here we can see that both c n 's and rf n 's are explicitly specified and the solution to the 
inhomogeneous GS equation associated with the multipolar boundary conditions is uniquely 
determined. The results in this Appendix establish the basis for further investigations of the 
evolution of the whole system in the main text. 



-24- 



REFERENCES 

Abdo, A. A., et al. 2011, Science, 331, 739 

Aly, J. J., 1984, ApJ, 283, 349 

Aly, J. J., 1991, ApJL, 375, 61 

Beloborodov, A. M. 2009, ApJ, 703, 1044 

Chen, J., et al. 1997, ApJL, 490, 191 

Dere, K. P., et al. 1999, ApJ, 516, 415 

Duncan, R. C, & Thompson, C, 1992, ApJL, 392, 9 

Feroci, M., Hurley, K., Duncan, R. C, & Thompson, C. 2001, ApJ, 549, 1021 

Forbes, T. G. Heliophysics: Space Storms and Radiation: Causes and Effects, 2010, Edited 
by Carolus, J. S. & George, L. S. (Cambridge Univ. Press), 159 

Forbes, T. G., & Isenberg, P. A. 1991, ApJ, 373, 294 

Forbes, T. G., & Priest, E. R. 1995, ApJ, 446, 377 

Gaensler, B. JVL, et al. 2005, Nature, 434,104 

Gavriil, F. P., Kaspi, V. M., Woods, P. M. 2002, Nature, 419, 142 

Gill, R., & Heyl, J. S. 2010, MNRAS, 407, 1926 

Gotz, D., Mereghetti, S., & Hurley, K. 2007, Ap&SS, 308, 51 

Isenberg, P. A., Forbes, T. G., & Demoulin, P. 1993, ApJ, 417, 368 

Jones, P. B. 2003, ApJ, 595, 342 

Kluzniak, W. & Ruderman, M., 1998, ApJL, 505, 113-117 
Komissarov, S. S., Barkov, M. & Lyutikov M., 2007, MNRAS, 374, 415 
Kouveliotou, C, et al. 1998, Nature, 393, 235 
Levin & Lyutikov, larXiv: 1204.26051 

Lin, J., Forbes, T. C, Isenberg, P. A., & Demoulin, P. 1998, ApJ, 504, 1006 



-25 - 

Lundquist, S. 1950, Ark. Fys., 2, 361 
Lyutikov, M. 2006, MNRAS, 367, 1602 
Mazets, E. P., et al. 1979, Nature, 282, 587 
Mereghetti, S. 2008, A&AR, 15, 225 
Mereghetti, S. k Stella, L. 1995, ApJL, 442, 17 
Palmer, D. M., et al. 2005, Nature, 434, 1107 

Pavan, L., Turolla, R., Zane, S., k Nobili, L. 2009, MNRAS 395, 753 
Perna, R., k Gotthelf, E. V. 2008, ApJ, 681, 522 

Priest E., k Forbes T. 2000, Magnetic Reconnection. MHD Theory and Applications. Cam- 
bridge Univ. Press, Cambridge 

Press, W. H., Teukolsky, S. A., Vetterling, W. T., k Flannery, B. P. 1992, Numerical Recipes 
in FORTRAN 2nd Ed., Cambrdige Univ. Press, Cambridge 

Ruderman, M. 1991, ApJ, 366, 261 

Shafranov, V. D. 1966, Rev. Plasma Phys., 2, 103 

Sturrock, P. A. 1991, ApJ, 380, 655 

Thompson, C, k Duncan, R. C. 1996, ApJ, 473, 322 

Thompson, C, k Duncan, R. C. 2001, ApJ, 561, 980 

Thompson, C, Lyutikov, M., k Kulkarni, S. R. 2002, ApJ, 574, 332 

Uzdensky, D. A., Cerutti, B. k Begelman, M. C. 2011, ApJL, 737, 40 

van Ballegooijen, A. A., & Martens, P. C. H. 1989, ApJ, 343, 971 

Woods, P. M., et al., 2001, ApJ, 552, 748 

Woods, P. M., k Thompson, C, 2006, in Compact Stellar X-Ray Sources, ed. W. H. G. 
Lewin k van der Klis (Cambridge Univ. Press), 547 

Woods, P. M., et al., 2007, ApJ, 654, 470 

Yu, C, 2011a, MNRAS, 411, 2461 



-26 - 



Yu, C, 2011b, ApJ, 738, 75 



This preprint was prepared with the A AS IATgX macros v5.2. 



